#include "projector.h"
#include "symconst.h"
#include "misc.h"
#include "utils.h"
using namespace cnst;

template <typename T>
bool find_edge(T edis,Graph& Gtmp)
{
  Graph::edge_iterator ei,ei_end;
  for(tie(ei,ei_end) = edges(Gtmp);ei != ei_end;++ei)
    if(*ei == Gtmp.global_to_local(edis))return true;
  return false;
}


void graph_ex_subs(exmap subs_to,Graph& G)
{
  Graph::edge_iterator ei,ei_end;
  Graph::vertex_iterator vi,vi_end;
   // substituting  new expressions for external momentum in edges and vertices
  for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
    {
      cout<<G[*ei].gmomentum<<"  "<<G[*ei].gmomentum.subs(subs_to)<<endl;
      cout<<"expr "<<G[*ei].gexpr<<"   "<<G[*ei].gexpr.subs(subs_to)<<endl;
      G[*ei].gmomentum = G[*ei].gmomentum.subs(subs_to);
      G[*ei].gexpr = G[*ei].gexpr.subs(subs_to);
   }
  ///vertices
  for(tie(vi,vi_end) = vertices(G);vi!=vi_end;++vi)
    {
      cout<<"expr "<<G[*vi].gexpr<<"   "<<G[*vi].gexpr.subs(subs_to)<<endl;
      G[*vi].gexpr = G[*vi].gexpr.subs(subs_to);
    }
}

void local_mom_subs(Graph& G,dia_ext_legs e_legs)
{
  exmap sm;
  sm[i1]=0;
  sm[i2]=0;
  sm[k1]=0;
  sm[k2]=0;
  graph_ex_subs(sm,G);
}

/// get summ of momentums external localy or globaly to VDIS,vdis is local descriptor
ex get_ext_mom_of_vertex(Graph::vertex_descriptor vdis,dia_ext_legs e_legs,Graph& Gtmp)
{
  Graph::vertex_iterator viu,viu_end,vi,vi_end;
  Graph::out_edge_iterator oei,oei_end;
  Graph::edge_descriptor edis;
  ex ret_mom = 0; 
  if(e_legs.count(Gtmp.local_to_global(vdis))>0)
    ret_mom += e_legs[Gtmp.local_to_global(vdis)].front().first;
  cout<<"from global vertex: "<<ret_mom<<endl;
  Graph G = Gtmp.root();
  ex fp = 0;
  for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(vdis),G);oei != oei_end;++oei)
    {
      if(G[*oei].from_vertex == Gtmp.local_to_global(vdis)) 
        {
          cout<<"out"<<G[*oei].gmomentum<<endl;
          fp-=G[*oei].gmomentum;
        }
      else if(G[*oei].to_vertex == Gtmp.local_to_global(vdis)) 
        {
          cout<<"in"<<G[*oei].gmomentum<<endl;
          fp+=G[*oei].gmomentum;
        }
    }
  cout<<"glob mom scaned "<<fp<<endl;
  ex fpl =0;
  for(tie(oei,oei_end) = out_edges(vertex(vdis,Gtmp),Gtmp);oei!=oei_end;++oei)
    {
      if(Gtmp[*oei].from_vertex ==Gtmp.local_to_global(vdis)) fpl+=Gtmp[*oei].gmomentum;
      else if(Gtmp[*oei].to_vertex==Gtmp.local_to_global(vdis)) fpl-=Gtmp[*oei].gmomentum;
    }
  cout<<"Local mom: "<<fpl<<endl;
  ret_mom+=(fp - fpl);
  cout<<"  MOMENTUM: "<<ret_mom<<endl;
  return ret_mom;
}


exvector get_input_momentums_vector(dia_ext_legs e_legs,Graph& G)
{
  exvector v1;
  Graph::vertex_iterator vi,vi_end;
  for(tie(vi,vi_end) = vertices(G);vi!=vi_end;++vi)
    {
      ex tn =  get_ext_mom_of_vertex(*vi,e_legs,G);
      if(tn!=0)  v1.push_back(tn);
      cout<<"VM: "<<tn<<endl;
    }
  return v1;
}

/// analog of Graph::find_vertex()
// synpsys:
//  edis to find,global, local graph
/*bool find_edge(Graph::edge_descriptor edis,Graph& Gtmp)
{
  cout<<" edge to find: "<<edis<<endl;
  Graph::edge_iterator ei,ei_end;
  for(tie(ei,ei_end) = edges(Gtmp);ei != ei_end;++ei)
    {
      cout<<"edges in graph "<<Gtmp.local_to_global(*ei)<<endl;
      if(Gtmp.local_to_global(*ei) == edis)return true;
     
    }
  return false;
}
*/

 ///!!!!!!!!!    HARD PROGRAMMING !!!!!!!!!!!!
 ///  find_edge implementation
bool find_edge(Graph::edge_descriptor edis, Graph& G)
{
  Graph::vertex_iterator vi,vi_end;
  Graph::edge_iterator ei,ei_end;
  cout<<"edge blyyaattt "<<edis<<"      "<<source(edis,G.root())<<"  "<<target(G.global_to_local(edis),G)<<endl;
  cout<<G.find_edge(edis).second<<endl;
  //  for(tie(vi,vi_end)=vertices(G);vi!=vi_end;++vi)

}

std::vector<field> get_local_ext_fields(Graph& Gtmp,dia_ext_legs e_l)
{
  print_graph(Gtmp);
  std::vector<field> field_vec;
  Graph G=Gtmp.root();
  print_graph(G);
  Graph::vertex_iterator vi,vi_end;
  Graph::edge_iterator eit,eit_end;
  //  for(tie(eit,eit_end)=edges(G);eit!=eit_end;++eit)
  // find_edge(*eit,Gtmp);

  for(  tie(vi,vi_end)=vertices(Gtmp);vi!=vi_end;++vi)
    {
      if( e_l.count( Gtmp.local_to_global(*vi))>0)
        {
          BOOST_FOREACH(mom_field_pair mfp_1,e_l[Gtmp.local_to_global(*vi)])
            {
              field_vec.push_back(mfp_1.second);
              cout<<"ext field " <<mfp_1.second<<endl;
            }
        }
      cout<<"Vertex "<<*vi<<" found "<<Gtmp.find_vertex(*vi).second<<endl;
      //      if((Gtmp.find_vertex(*vi).second)&&(out_degree(*vi,G)>out_degree(Gtmp.global_to_local(*vi),Gtmp)))
      if(out_degree(*vi,Gtmp)<out_degree(Gtmp.local_to_global(*vi),G))
        {
          Graph::out_edge_iterator oei,oei_end;
          Graph::edge_descriptor edis;
          for(tie(oei,oei_end)=out_edges(Gtmp.local_to_global(*vi),G);oei!=oei_end;++oei)
            //      if(!find_edge(*oei,Gtmp))
            //            cout<<"src: "<<source(*oei,G)<<" targ: "<<target(*oei,G)<<endl;
            if(!Gtmp.find_edge(*oei).second)
              {
                cout<<"edge "<<(*oei)<<" not found"<<endl;
                edis = *oei;
                if(G[edis].from_vertex==Gtmp.local_to_global(*vi))field_vec.push_back(G[*oei].from_field);
                if(G[edis].to_vertex==Gtmp.local_to_global(*vi))field_vec.push_back(G[*oei].to_field);
              }          
            else cout<<"edge in subgraph "<<(*oei)<<endl;
        }
    }
  sort(field_vec.begin(),field_vec.end());
  return field_vec;
}





tuple<exmap,exmap,exvector> get_local_mom_subs_map(exvector ext_vec)
///please mark all momentums before input,otherwise new moments are input only!!!
//ex afdcpp::subs_ext_mom(ex in_ex,exvector ext_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | External momentum substitution      | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  //symbol l1("l1"),l2("l2"),l3("l3"),l4("l4"),l5("l5");
  //const ex lm[] = {l1,l2,l3,l4,l5};
  int n = 0;
  //std::list<ex> m_lst_lst;
  //list of all momentum symbols in external lines 
  std::set<symbol,ex_is_less> m_set;
  BOOST_FOREACH(ex e478,ext_vec)
    {

      cout<<"mom of vert: "<<endl;
      std::copy(e478.preorder_begin(),e478.preorder_end(),std::ostream_iterator<ex>(cout, "\n"));
      cout<<"my way:"<<endl;
      
      for (const_preorder_iterator it = e478.preorder_begin();it!=e478.preorder_end();++it)
        {
          if(is_a<symbol>(*it))m_set.insert(ex_to<symbol>(*it));
          cout<<"---> "<<(*it)<<" is a symbol: "<<is_a<symbol>(*it)<<endl;
        }
        
          //      cout<<"solution subs "<<(*it)<<" = "<<lsolve(e478==lm[n],(*it))<<endl;
    }
  cout<<"not unique symbols: "<<m_set.size()<<endl;
  /*  BOOST_FOREACH(symbol e634,m_set)
    {
      cout<<"not unique "<<e634<<endl;
      }*/
  //  m_set.sort();
  // m_set.unique();
  BOOST_FOREACH(symbol e634,m_set)
    {
      cout<<"unique "<<e634<<endl;
    }
  /// creation of momentum index map
  //std::map<symbol,int> mom_idx;
  exhashmap<int> mom_idx;
  BOOST_FOREACH(ex e789,m_set)
    {
      BOOST_FOREACH(ex e45,ext_vec)
        {
          if(e45.has(e789))
            {
              cout<<e45<<" has "<<e789<<endl;
              mom_idx[ex_to<symbol>(e789)]++;
            }
        }
    }
  for(exhashmap<int>::iterator it = mom_idx.begin();it != mom_idx.end();++it)
    cout<<"mom: "<<it->first<<" mom_idx: "<<it->second<<endl;
  /// finding of minimal momentum index in each ext_mom
  /// and mom map creation
  //list in form (ext_mom,min_symb_mom,idx),where min is a momentum with minimal index
  std::vector< tuple<ex,symbol,int> > em_m_lst;
  BOOST_FOREACH(ex e45,ext_vec)
    {
      symbol buf_mom;
      int buf_idx = 1000;
      BOOST_FOREACH(symbol e789,m_set)
        {
          if(e45.has(e789) && (mom_idx[e789] < buf_idx))
            {
              buf_idx = mom_idx[e789];
              buf_mom = e789;
            }
        }
      em_m_lst.push_back(make_tuple(e45,buf_mom,buf_idx));
    }
  sort(em_m_lst.begin(),em_m_lst.end(),nops_less_comp);
  sort(em_m_lst.begin(),em_m_lst.end(),idx_less_comp);

  typedef tuple<ex,symbol,int> mom_tuple;
  
  exmap mom_subs_map,mom_subs_back_map;    
  exvector l_mom_vec;
  enum{ext_mom_expression,mom_symbol,mom_index};
  ex buf_momentum = 0;
  BOOST_FOREACH(mom_tuple m_pair,em_m_lst)
    {
      if(n<(em_m_lst.size()-1))
        {
          mom_subs_map[m_pair.get<mom_symbol>()] = lsolve(m_pair.get<0>()==lm[n],m_pair.get<mom_symbol>());
          mom_subs_back_map[lm[n]] = m_pair.get<ext_mom_expression>();
          cout<<"Ext mom: "<<m_pair.get<0>()<<" min "<<m_pair.get<1>()<< " index "<<mom_idx[m_pair.get<1>()]<<endl;
          l_mom_vec.push_back(lm[n]);
          buf_momentum+=lm[n];
        }
      else
        {
          mom_subs_map[m_pair.get<mom_symbol>()] = lsolve(m_pair.get<0>()==-buf_momentum,m_pair.get<mom_symbol>());
          //          mom_subs_back_map[lm[n]] = m_pair.get<ext_mom_expression>();
          cout<<"Last mom: "<<m_pair.get<0>()<<" min "<<m_pair.get<1>()<< " index "<<mom_idx[m_pair.get<1>()]<<endl;
          l_mom_vec.push_back(-buf_momentum);
        }
      n++;
    }
  //  cout<<colors::blue<<"SMOM:  "<<in_ex.subs(mom_subs_map)<<colors::def_col<<endl;
  return make_tuple(mom_subs_map,mom_subs_back_map,l_mom_vec);
}


/// Try use copy_graph from boost/copy.hpp
/// NEED REVIEW!!!! some fields may be empty!!!!!!!!!!
Graph copy_graph(Graph& G)
{
      Graph::vertex_iterator viu,viu_end;
      Graph::edge_iterator ei,ei_end;
      Graph::out_edge_iterator oei,oei_end;
      Graph GT;              /// Graph after copy
      Graph::edge_descriptor edis;
      Graph::vertex_descriptor gvdis;
      bool inserted;
      get_property(GT,graph_num) = get_property(G,graph_num);
      for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
        {
          tie(edis, inserted) = add_edge(G.local_to_global(source(*ei,G)),G.local_to_global(target(*ei,G)),GT);
              GT[edis].type = G[*ei].type;
              GT[edis].gmomentum = G[*ei].gmomentum;
              GT[edis].from_vertex = G[*ei].from_vertex;
              GT[edis].to_vertex = G[*ei].to_vertex;
              GT[edis].gexpr = G[*ei].gexpr;
              cout<<"Inserted full edge"<<endl;
        }
      for (tie(viu, viu_end) = vertices(GT);viu != viu_end; ++viu)
        {
          gvdis = vertex(*viu,G); // vertex discriptor of root graph vert
          //  gvdis = G.global_to_local(*viu);
          GT[*viu].gexpr = G.root()[gvdis].gexpr ;
          GT[*viu].gmomentum = G.root()[gvdis].gmomentum ;
          GT[*viu].external = G.root()[gvdis].external;
          cout<<"GT has "<<num_vertices(GT)<<" vertices\n";
        
        }
      return GT;
}




/*Graph get_clear_graph(Graph& Gtmp,dia_ext_legs e_legs)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |      Constructing T(G) terms        | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  if(Gtmp.is_root())
    cout<<"\n ERROR: Can't apply T(G) on G, where G is a root graph object "<<endl;
  else
    {
      cout<<"and Gtmp"<<endl;
      print(Gtmp);
      Graph::vertex_descriptor gvdis,vdis;
      bool is_in_graph;
      Graph::vertex_iterator viu,viu_end;
      Graph::edge_iterator ei,ei_end;
      Graph::out_edge_iterator oei,oei_end;
      Graph G = Gtmp.root(); /// Main diagram graph

      cout<<"print Gtmp.root"<<endl;
      print(G);

      Graph GT;              /// Graph after T-operation
      get_property(GT,graph_num) = get_property(Gtmp,graph_num);
      bool from_is_in_div_graph,to_is_in_div_graph;
      bool inserted;
      Graph::edge_descriptor edis;
      // tuple<vdis,edis,field,momentum>
      tuple_lst tleg_lst;
      for (tie(ei, ei_end) = edges(G);ei != ei_end; ++ei)
        {
          //          tie(vdis,from_is_in_div_graph) = Gtmp.find_vertex(G[*ei].from_vertex);
          //          tie(vdis,to_is_in_div_graph) = Gtmp.find_vertex(G[*ei].to_vertex);

          from_is_in_div_graph =find_vert(G[*ei].from_vertex,Gtmp);
          to_is_in_div_graph = find_vert(G[*ei].to_vertex,Gtmp);
          
          if(from_is_in_div_graph && to_is_in_div_graph)                /// edge is in div subgraph
            {
                            cout<<"full edge in subgraph"<<endl;
            }
          if((!from_is_in_div_graph) && (!to_is_in_div_graph))            /// adding edge to T(G)
            {
              tie(edis, inserted) = add_edge(G[*ei].from_vertex,G[*ei].to_vertex,GT);
              GT[edis].type = G[*ei].type;
              GT[edis].gmomentum = G[*ei].gmomentum;
              GT[edis].from_vertex = G[*ei].from_vertex;
              GT[edis].to_vertex = G[*ei].to_vertex;
              GT[edis].gexpr = G[*ei].gexpr;
              cout<<"Inserted full edge"<<endl;
              
            }
          if((!from_is_in_div_graph) && to_is_in_div_graph)           /// from field     
            {
              tleg_lst.push_back(make_tuple(G[*ei].from_vertex,*ei,G[*ei].from_field,G[*ei].gmomentum));
              cout<<"To "<<G[*ei].to_vertex<<"is in div graph, and adding"<<G[*ei].from_vertex<<endl;
            }
          if(from_is_in_div_graph && (!to_is_in_div_graph))         /// to field     
            {
              tleg_lst.push_back(make_tuple(G[*ei].to_vertex,*ei,G[*ei].to_field,G[*ei].gmomentum));
              cout<<"From "<<G[*ei].from_vertex<<"is in div graph, and adding"<<G[*ei].to_vertex<<endl;
            }
        }
      /// Adding vertex information
      /// Check if no vertices in GT
      ///
      /// !!!!!!!!!!SINGLE VERTICES ARE NOT ADDED!!!!!!!!!!!!
      ///
      /// Adding info in all GT vertices
      for (tie(viu, viu_end) = vertices(GT);viu != viu_end; ++viu)
        {
          gvdis = vertex(*viu,G); // vertex discriptor of root graph vert
          GT[gvdis].gexpr = G[gvdis].gexpr ;
          GT[gvdis].gmomentum = G[gvdis].gmomentum ;
          GT[gvdis].external = G[gvdis].external;
          cout<<"GT has "<<num_vertices(GT)<<" vertices\n";
          
        }
      ext_vert_mom_field_tuple_lst l3 = get_ext_legs(Gtmp);
      cout<<"FNumber of added verts: "<<tleg_lst.size()<<" and ext in GT "<<num_ext_vert(GT)<<" and real ext in M = "<<l3.size()<<" verts to shrink: "<< l3.size()+tleg_lst.size()<<endl;
      /// number of external legs
      int n_ext_l = l3.size();
      /// number of internal legs
      int n_int_l = tleg_lst.size(); 
      /// number of legs in subgraph to be shrinked
      int n_legs = n_ext_l+n_int_l;
      switch(n_legs)
        {
        case 2: // ins propagator
          {
            cout<<"\n\n  ext legs: "<<n_ext_l<<"       I want propagator!!!\n\n"<<endl;
            
          }
        case 3: // ins 3 part vertex
          {
             cout<<"\n\n  ext legs: "<<n_ext_l<<"       I want 3 part vert!!!\n\n"<<endl;
          }
        case 4: // ins 4 part vertex
          {
            vdis = Gtmp.local_to_global(vertex(0,Gtmp));
            cout<<" VDIS:"<<vdis<<endl;
            if(n_int_l > 0)   /// if edges between new vertex and old exists
              for(tuple_lst::iterator titer = tleg_lst.begin();
                  titer != tleg_lst.end();++titer)
                {
                  cout<<"Before ins"<<endl;
                  print_edges2(GT, get(vertex_index, GT), get(edge_index, GT));
                  print_vertices(GT,get(vertex_index, GT));
                  //print_edges(GT,get(edge_index, GT));
                  print_graph(GT);
                  /// add edge from new vertex (vdis) to vertex <0> in half of edge
                  tie(edis, inserted) = add_edge(vdis,get<0>(*titer),GT);
                  if(get<0>(*titer) == G[get<1>(*titer)].from_vertex)//edge from M(G)
                    {
                      tie(GT[edis].type,GT[edis].gmomentum,GT[edis].gexpr)=
                        T_expr_prop(G[get<1>(*titer)].from_field,
                                    G[get<1>(*titer)].to_field,
                                    G[get<1>(*titer)].to_vertex,vdis,
                                    G[get<1>(*titer)].gmomentum);
                      GT[edis].from_vertex = get<0>(*titer);
                      GT[edis].to_vertex = vdis;
                      GT[edis].from_field = G[get<1>(*titer)].from_field;
                      GT[edis].to_field = G[get<1>(*titer)].to_field;
                      /// adding expression to vertex from root graph
                      GT[get<0>(*titer)].gexpr = G[get<0>(*titer)].gexpr;
                      GT[get<0>(*titer)].external = G[get<0>(*titer)].external;
                      GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
                    }
                  if(get<0>(*titer) == G[get<1>(*titer)].to_vertex) // edge to M(G)
                    {
                      tie(GT[edis].type,GT[edis].gmomentum,GT[edis].gexpr)=
                        T_expr_prop(G[get<1>(*titer)].from_field,
                                    G[get<1>(*titer)].to_field, vdis,
                                    G[get<1>(*titer)].to_vertex,
                                    G[get<1>(*titer)].gmomentum);
                      GT[edis].from_vertex = vdis;
                      GT[edis].to_vertex = get<0>(*titer);
                      GT[edis].from_field = G[get<1>(*titer)].from_field;
                      GT[edis].to_field = G[get<1>(*titer)].to_field;
                      /// adding expression to vertex from root graph
                      GT[get<0>(*titer)].gexpr = G[get<0>(*titer)].gexpr;
                      GT[get<0>(*titer)].external = G[get<0>(*titer)].external;
                      GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
                    }
                  
                  cout<<"After ins"<<endl;
                  print_vertices(GT,get(vertex_index, GT));
                  //print_edges(GT,get(edge_index, GT));
                  print_graph(GT);
                  cout<<"-----"<<endl;
                }
            else        /// new edges are no added diagram==single vertex
              {
               vdis =  add_vertex(GT);
               cout<<"Single vert "<<endl;
                  print_graph(GT);
                  cout<<"-----"<<endl;
              }
            
            // if external legs exists in, than it is external vertex
            GT[vdis].external = (n_ext_l > 0);
            tie(GT[vdis].gmomentum,GT[vdis].gexpr)=T_expr4(vdis,tleg_lst,l3);
            //GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
            cout<<ex_from_graph(GT)<<endl;

            /*tie(edis, inserted) = add_edge(G[*ei].from_vertex,G[*ei].to_vertex,GT);
              GT[edis].type = G[*ei].type;
              GT[edis].gmomentum = G[*ei].gmomentum;
              GT[edis].from_vertex = G[*ei].from_vertex;
              GT[edis].to_vertex = G[*ei].to_vertex;
              GT[edis].gexpr = G[*ei].gexpr;*/
/*
          }
        }
      cout<<"print T(G)"<<endl;
      print(GT);
        
      cout<<"Creates T(G)"<<endl;
      return GT;
    }
      
            

/*
  //  cout<<"vertieces in parent :" <<num_vertices(g.parent())<<"And in child:"<<num_vertices(g)<<endl;
  /// Marking vertieces as external for subgraph
      for (tie(viu, viu_end) = vertices(Gtmp);viu != viu_end; ++viu) 
        {
          if(out_degree(*viu,Gtmp) < out_degree( Gtmp.local_to_global(*viu), G) )
            {
              Gtmp[vertex(*viu,Gtmp)].external = true;
               
              cout<<"################external found!!!!"<<endl;
              /// Vertex Momentum Flow Calculation
              ex fp;
              for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(*viu),G);oei != oei_end;++oei)
                fp+=G[*oei].gmomentum;
              for(tie(oei,oei_end) = out_edges(vertex(*viu,Gtmp),Gtmp);oei!=oei_end;++oei)
                fp-=Gtmp[*oei].gmomentum;
              cout<<"  MOMENTUM: "<<fp<<endl;
              */
/*               ex e37 = ex_from_graph(Gtmp);
  symbol n("n");
/// Working example
 cout<< "\nExternal changed 1 " << e37.subs(k1-p1+wild() == n + wild())<<endl;
 
 
  cout<< "\nExternal changed 2 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
   cout<< "\nExternal changed 3 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
  

/*
              Gtmp[*viu].gmomentum = fp;
            }
        }
  */
/*
  //ex mom = k1+p1;
  symbol k;
  ex mom = k;
  symbol sim = m0;

}

*/







ex projector(Graph& G,exmap back_map,std::vector<field> f_vec,exvector leg_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Expressioning numerator in terms    | ///
/// | of scalar formfactors,momentum subs | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  
  //      std::vector<field> lf = get_ext_fields(G);
  //    sort(lf.begin(),lf.end());
      BOOST_FOREACH(field f23,f_vec)
        {
          cout<<"ext field "<<f23<<endl;
        }
      switch(f_vec.size())
        {
        case 2:
          {
            if(f_vec[0]==e && f_vec[1]==E)cout<<"QED e self energy found"<<endl;
            if(f_vec[0]==A && f_vec[1]==A)cout<<"QED polarization found"<<endl;
            return 0;
          }
        case 3:
          {
            if(f_vec[0]==e && f_vec[1]==E && f_vec[2]==A)
              {
                cout<<"QED vert found"<<endl;
                /*               ex out_e = ex_from_graph(G);
                //exvector 
                std::vector<ex> indvector = out_e.get_free_indices();
                if(indvector.size() == 1)
                  {
                    cout<<indvector[0]<<endl;
                    varidx vi_in;
                    if(is_a<varidx>(indvector[0]))
                     vi_in = ex_to<varidx>(indvector[0]);
                    else cout<<"BAAAAAAD"<<endl;
                    out_e = dirac_trace(dirac_gamma(vi_in.toggle_variance())*out_e).simplify_indexed();
                    cout<<out_e<<endl;
                    exvector ev = get_ext_momentum(G);
                    exvector lum = loop_momentum(out_e,ev);
                
                    return out_e;*/
                  }
                else 
                  {
                    cout<<colors::red<<"ERROR: more then one index in QED vertex"<<colors::def_col<<endl;
                    return 0;
                  }
              }
        case 4:
          {
            if(f_vec[0]==phi && f_vec[1]==phi && f_vec[2]==phi && f_vec[3] == phi)
              {
                cout<<"phi4 vert found"<<endl;
                /*
                  ex out_e = ex_from_graph(G);
                  cout<<"Scalarized: "<<scalarize_vectors(out_e)<<endl;
                  exvector ev1 = get_ext_momentum(G);
                  exvector lum = loop_momentum(out_e,ev1);
                  ex soue = subs_ext_mom(out_e,ev1);
                  feynman(soue,lum,ev1);
                  /// GET_MASSES function!!!
                  /// from graph (data in propagator or find in expression)
                  //MSbar not depend on momentum and other
                  // One-loop integral function
                  symbol rho("rho");
                  ex C = Scalar2Pt(p1, mS, mS, 1, 1, rho);//, eps);
                  // Perform epsilon-expansion
                  C = C.series(eps == 0, 1);
                  // Extract the divergent and finite parts
                  ex div = out_e.numer()*1/eps*C.coeff(eps, -1).subs(rho == 0).normal();
                  
                  return div;*/
              }
          }
          /*        default:
                    {
                    cout<<" ERROR: Unknown number of external lines "<<lf.size()<<endl;
                    return 0;
                    }*/
        }
      symbol ret_sym;
      return ret_sym;
}

ex root_projector(Graph& G,std::vector<field> f_vec,exvector leg_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Expressioning numerator in terms    | ///
/// | of scalar formfactors,momentum subs | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  
  //      std::vector<field> lf = get_ext_fields(G);
  //    sort(lf.begin(),lf.end());
      BOOST_FOREACH(field f23,f_vec)
        {
          cout<<"ext field "<<f23<<endl;
        }
      switch(f_vec.size())
        {
        case 2:
          {
            if(f_vec[0]==e && f_vec[1]==E)cout<<"QED e self energy found"<<endl;
            if(f_vec[0]==A && f_vec[1]==A)cout<<"QED polarization found"<<endl;
            return 0;
          }
        case 3:
          {
            if(f_vec[0]==e && f_vec[1]==E && f_vec[2]==A)
              {
                cout<<"QED vert found"<<endl;
                /*               ex out_e = ex_from_graph(G);
                //exvector 
                std::vector<ex> indvector = out_e.get_free_indices();
                if(indvector.size() == 1)
                  {
                    cout<<indvector[0]<<endl;
                    varidx vi_in;
                    if(is_a<varidx>(indvector[0]))
                     vi_in = ex_to<varidx>(indvector[0]);
                    else cout<<"BAAAAAAD"<<endl;
                    out_e = dirac_trace(dirac_gamma(vi_in.toggle_variance())*out_e).simplify_indexed();
                    cout<<out_e<<endl;
                    exvector ev = get_ext_momentum(G);
                    exvector lum = loop_momentum(out_e,ev);
                
                    return out_e;*/
                  }
                else 
                  {
                    cout<<colors::red<<"ERROR: more then one index in QED vertex"<<colors::def_col<<endl;
                    return 0;
                  }
              }
        case 4:
          {
            if(f_vec[0]==phi && f_vec[1]==phi && f_vec[2]==phi && f_vec[3] == phi)
              {
                cout<<"phi4 vert found"<<endl;
                /*
                  ex out_e = ex_from_graph(G);
                  cout<<"Scalarized: "<<scalarize_vectors(out_e)<<endl;
                  exvector ev1 = get_ext_momentum(G);
                  exvector lum = loop_momentum(out_e,ev1);
                  ex soue = subs_ext_mom(out_e,ev1);
                  feynman(soue,lum,ev1);
                  /// GET_MASSES function!!!
                  /// from graph (data in propagator or find in expression)
                  //MSbar not depend on momentum and other
                  // One-loop integral function
                  symbol rho("rho");
                  ex C = Scalar2Pt(p1, mS, mS, 1, 1, rho);//, eps);
                  // Perform epsilon-expansion
                  C = C.series(eps == 0, 1);
                  // Extract the divergent and finite parts
                  ex div = out_e.numer()*1/eps*C.coeff(eps, -1).subs(rho == 0).normal();
                  
                  return div;*/
              }
          }
          /*        default:
                    {
                    cout<<" ERROR: Unknown number of external lines "<<lf.size()<<endl;
                    return 0;
                    }*/
        }
      symbol ret_sym;
      return ret_sym;
}


void get_shrinked_graph(Graph& G,Graph& subG,ex subs_ex)
{
  Graph out_G;
  Graph::edge_iterator ei,ei_end;
  //  print(G);
  // print(subG);
  // if edge is in subG -not added, if connects G and subG, then added frm G vertex to shrinked V
  Graph::vertex_descriptor vertex_shrink_to = subG.local_to_global(0); // first vertex of shrinked subgraph
  cout<<"debug"<<endl;
  Graph::edge_descriptor edis;
  bool inserted;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    {
      cout<<"edge "<<*ei<<" found= "<<subG.find_edge(G.local_to_global(*ei)).second<<endl;
      if(!subG.find_edge(G.local_to_global(*ei)).second) // edges not in subG
        {
          if(subG.find_vertex(G.local_to_global(source(*ei,G))).second) // edge source is in subG
            {
              tie(edis,inserted) = add_edge(vertex_shrink_to,G.local_to_global(target(*ei,G)),out_G);
              out_G[edis].from_vertex = vertex_shrink_to; ///need review
              out_G[edis].to_vertex =G.local_to_global(target(*ei,G));//need review
              out_G[G.local_to_global(target(*ei,G))].external = G[target(*ei,G)].external;
              out_G[G.local_to_global(target(*ei,G))].gexpr = G[target(*ei,G)].gexpr;
            }
          else if(subG.find_vertex(G.local_to_global(target(*ei,G))).second) // edge target is in subG
            {
              tie(edis,inserted) = add_edge(G.local_to_global(source(*ei,G)),vertex_shrink_to,out_G);
              out_G[edis].from_vertex = G.local_to_global(source(*ei,G));//need review
              out_G[edis].to_vertex =vertex_shrink_to; ///need review
              out_G[G.local_to_global(source(*ei,G))].external = G[source(*ei,G)].external;
              out_G[G.local_to_global(source(*ei,G))].gexpr = G[source(*ei,G)].gexpr;
            }
          else // not u,no v from (u,v) are in subG
            {
              tie(edis,inserted) = add_edge(G.local_to_global(source(*ei,G)),G.local_to_global(target(*ei,G)),out_G);
              out_G[edis].from_vertex = G.local_to_global(source(*ei,G));//need review
              out_G[edis].to_vertex = G.local_to_global(target(*ei,G));//need review
              out_G[G.local_to_global(source(*ei,G))].external = G[source(*ei,G)].external;
              out_G[G.local_to_global(target(*ei,G))].external = G[target(*ei,G)].external;
              out_G[G.local_to_global(source(*ei,G))].gexpr = G[source(*ei,G)].gexpr;
              out_G[G.local_to_global(target(*ei,G))].gexpr = G[target(*ei,G)].gexpr;
            }
          out_G[edis].type = G[*ei].type;
          out_G[edis].gmomentum = G[*ei].gmomentum;
          out_G[edis].gexpr = G[*ei].gexpr;
          out_G[edis].from_field = G[*ei].from_field;
          out_G[edis].to_field = G[*ei].to_field;
        }
    }
  out_G[vertex_shrink_to].external = false;
  out_G[vertex_shrink_to].gexpr = subs_ex;
  print(out_G);
}
